In this analysis, we use R to read in LandSat8 data, create a vegitation index and calculate zonal statistics of a raster. The data come from the 2006 National Land Cover Database.

First we relcassify the NLCD data into two classes, based on the value of the raster. In this case, I want to classify the pixels as to whether they are developed or undeveloped.

We then use a census tract polygon layer to calculate the proportion of each tract’s land area that is developed vs undeveloped.

Raster read in

library(rgdal)
## Loading required package: sp
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/ozd504/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/ozd504/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/ozd504/R/win-library/4.0/rgdal/proj
library(raster)
library(gdalUtils)
nlcd<-raster("~/OneDrive - University of Texas at San Antonio/classes/dem7093/dem7093_21//data/nlcd16_48_lc/nlcd_2016_land_cover_48_20190424.img.vat.img")

 nlcdproj<-gdalwarp(srcfile ="~/OneDrive - University of Texas at San Antonio/classes/gis_classwork/nlcd16_48_lc/nlcd_2016_land_cover_change_index_48_20190424.img.vat.img", dstfile = "~/OneDrive - University of Texas at San Antonio/classes/gis_classwork/outputnlcd.tif",t_srs = nlcd , output_Raster = T)
## Warning in system(cmd, intern = TRUE): running command '"C:
## \OSGeo4W64\bin\gdalwarp.exe" -t_srs "raster(ncol=44253, nrow=41260,
## xmn=-1075365, xmx=252225, ymn=309855, ymx=1547655, crs='+proj=aea
## +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0
## +datum=WGS84 +units=m +no_defs')" -of "GTiff" "~/OneDrive -
## University of Texas at San Antonio/classes/gis_classwork/nlcd16_48_lc/
## nlcd_2016_land_cover_change_index_48_20190424.img.vat.img" "~/OneDrive -
## University of Texas at San Antonio/classes/gis_classwork/outputnlcd.tif"' had
## status 1
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition
 raster::plot(nlcdproj, main="Original Raster image")

Clip raster to polygon

library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
## 
## Attaching package: 'sf'
## The following object is masked from 'package:gdalUtils':
## 
##     gdal_rasterize
library(tigris)
## To enable 
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
library(terra)
## terra version 1.1.4
## 
## Attaching package: 'terra'
## The following object is masked from 'package:rgdal':
## 
##     project
options(tigris_class = "sf")
counties<-tracts(state="TX", county = c("Bexar", "Comal"), cb=T, year=2016)
counties<-st_transform(counties, crs=2278)

# doesn't work

est<-terra::crop(nlcdproj, counties)
plot(est)

est2<-terra::mask(est, counties)
plot(est2)

# est<-mask(x=nlcdproj,mask=as_Spatial(counties))
# est2<-crop(x=nlcdproj, y=extent(as_Spatial(counties)))
# plot(est2)

Raster reclassification:

m<-c(1, 12, 0,13, 24, 1, 25, 95, 0)
rclmat<-matrix(m, ncol=3, byrow = T)
nlcdreclass<-raster::reclassify(x = est2, rcl = rclmat)

plot(nlcdreclass, main="Reclassified Raster image \n 1=developed, 0 = undeveloped",col= gray(100:0 / 100))

Load the census tract polygons

library(tigris)
satract<-tracts(state="48", county="029", year=2010)

library(sf)
satract<-st_as_sf(satract)
satract<-st_transform(satract, crs=  2278)

satract<-as(satract, "Spatial")
summary(satract)

Calculate zonal statistics

Now we find the mean of the binary reclassified data, to estimate the percent of each tract’s area that is developed.

First we reproject the raster to the same coordinate system as the census data. I get the projection string from ^^^ above.

nlcdreclass<-projectRaster(nlcdreclass, crs="+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")

satract$pct_developed<-as.numeric(extract(nlcdreclass, satract, fun=mean)) #this takes a minute or two
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster
satract<-st_as_sf(satract)

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     collapse, desc, intersect, near, select, union
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
satract<-satract%>%
  filter(complete.cases(pct_developed))

tractlines<-st_boundary(satract)
library(tmap)
tm_shape(satract)+
  tm_polygons("pct_developed",style="quantile", n=5,legend.hist=T )+
  tm_shape(tractlines)+
  tm_lines(col="black")+
  tm_format("World",
            title="San Antonio Developed Land Area",
            legend.outside=TRUE )

## [1] "#FFF8C4" "#FEDE86" "#FEAA38" "#EA6E13" "#B74202"
library(classInt)
m1<-satract%>%
  mutate(devcut=cut(satract$pct_developed,breaks=data.frame(classIntervals(var=satract$pct_developed, n=5, style="jenks")[2])[,1], include.lowest = T))%>%
  ggplot(aes( fill=devcut))+geom_sf()+
  scale_fill_brewer(palette = "Blues") + 
  scale_color_brewer(palette = "Blues")+
  geom_sf(data=tractlines,fill=NA, color="black")+
  guides(fill=guide_legend(title="Percent Developed"))+ggtitle(label="Percent of Land Area Developed, 2006", subtitle ="NLCD 2006")+theme(axis.text.x = element_blank(), axis.text.y = element_blank())+theme( legend.text = element_text(size = 14), plot.title = element_text(size = 18), strip.text.x = element_text(size = 16))

m1

I used this method in a paper I wrote in 2013, although in that paper I did the diversity of the land cover. That could be done like:

#nlcdproj<-projectRaster(nlcd, crs="+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")

library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
shannonVegan <- function(x, ...) {
    diversity(table(x), index="shannon")
}

#shanVegOut <- focal(nlcdproj, fun=shannon, pad=T)  

satract$var_land<-as.numeric(extract(nlcdproj, satract, fun=shannonVegan)) #this takes a minute or two
shannon<-tm_shape(satract)+
  tm_polygons("var_land",style="quantile", n=5,legend.hist=T )+
  tm_shape(tractlines)+
  tm_lines(col="black")+
  tm_format("World",
            title="San Antonio Shannon Diversity of Land Area",
            legend.outside=TRUE )

shannon

## [1] "#FFF8C4" "#FEDE86" "#FEAA38" "#EA6E13" "#B74202"

Exporting images for posters/presentations

I recommend using the ggsave() function to export images made by ggplot for use in other documents. For tmap, you can use tmap_save()

#ggsave(m1, path = "", filename = "sa_pct_devel.png", units = "in", width=8, height = 10)
tmap_save(shannon,filename = "~/OneDrive - University of Texas at San Antonio/classes/gis_classwork/shannon_tract.png", width = 8, height = 10, units = "in")

And you can insert that image into a doc or presentation.

Landsat imagery

Good simple description of landsat here

library(raster)
library(rgdal)
library(rgeos)
## rgeos version: 0.5-5, (SVN revision 640)
##  GEOS runtime version: 3.8.0-CAPI-1.13.1 
##  Linking to sp version: 1.4-2 
##  Polygon checking: TRUE
library(RColorBrewer)


# turn off factors
options(stringsAsFactors = FALSE)
landsatstack<- list.files("~/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/", 
  pattern=glob2rx("*B*.TIF$"), 
  full.names = T)

landsatstack
##  [1] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B1.TIF" 
##  [2] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B10.TIF"
##  [3] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B11.TIF"
##  [4] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B2.TIF" 
##  [5] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B3.TIF" 
##  [6] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B4.TIF" 
##  [7] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B5.TIF" 
##  [8] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B6.TIF" 
##  [9] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B7.TIF" 
## [10] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B8.TIF" 
## [11] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B9.TIF" 
## [12] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/classes/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_BQA.TIF"
lsb2<-raster(landsatstack[2])

plot(lsb2,col = gray(0:100 / 100))  

  rm(lsb2)      ; gc()
##            used  (Mb) gc trigger  (Mb)  max used   (Mb)
## Ncells  4441992 237.3    8505014 454.3   8505014  454.3
## Vcells 26776359 204.3  102689378 783.5 458129635 3495.3

Make raster stack of bands 1 - 7

lsstack<-stack(landsatstack[c(1,4:9)])

lsbrick<-brick(lsstack)
names(lsbrick)
## [1] "LC08_L1TP_027040_20200311_20200325_01_T1_B1"
## [2] "LC08_L1TP_027040_20200311_20200325_01_T1_B2"
## [3] "LC08_L1TP_027040_20200311_20200325_01_T1_B3"
## [4] "LC08_L1TP_027040_20200311_20200325_01_T1_B4"
## [5] "LC08_L1TP_027040_20200311_20200325_01_T1_B5"
## [6] "LC08_L1TP_027040_20200311_20200325_01_T1_B6"
## [7] "LC08_L1TP_027040_20200311_20200325_01_T1_B7"
names(lsbrick) <- gsub(pattern = "LC08_L1TP_027040_20200311_20200325_01_T1_",
                       replacement = "",
                       names(lsbrick))
names(lsbrick)
## [1] "B1" "B2" "B3" "B4" "B5" "B6" "B7"
ndvi<-(lsbrick[[5]] - lsbrick[[4]]) / (lsbrick[[5]] + lsbrick[[4]])
plot(ndvi,
     main = "NDVI around San Antonio",
     axes = FALSE, box = FALSE)       

hist(ndvi)

library(mapview)
mapview(nlcdreclass)
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 9494355 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  9494355 '
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
writeRaster(x= ndvi, 
    filename = "~/OneDrive - University of Texas at San Antonio/classes/gis_classwork/ndvi1.tif", format="GTiff", datatype="INT2S", overwrite=T)

Landsat plots

thermavg<-
plotRGB(lsbrick,
r = 6, g = 5, b = 2,
stretch = "lin",
axes = TRUE,
main = "RGB composite image\n showing Agricultural land - Landsat Bands 6, 5, 2")
box(col = "white")

plotRGB(lsbrick,
r = 5, g = 4, b = 3,
stretch = "lin",
axes = TRUE,
main = "Color infrared composite image\n Landsat Bands 5, 4, 3")
box(col = "white")